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Many biological and soft matter processes occur at high speeds in complex 3D environments, and 
developing imaging techniques capable of elucidating their dynamics is an outstanding experimental 
challenge. Here, we introduce Fourier Synthesis Optical Diffraction Tomography (FS-ODT), a novel 
approach for high-speed quantitative phase imaging capable of recording the 3D refractive index 
at kilohertz rates. FS-ODT introduces new pattern generation and inverse computational strate- 
gies that multiplex tens of illumination angles in a single tomogram, dramatically increasing the 
volumetric imaging rate. We validate FS-ODT performance by imaging samples of known compo- 
sition and accurately recovering the refractive index for increasing pattern complexity. We further 
demonstrate the capabilities of FS-ODT for probing complex systems by studying the hindered dif- 
fusion of colloids in solution and the motility of single-cellular bacterial swimmers. We believe that 
FS-ODT is a promising approach for unlocking challenging imaging regimes in biophysics and soft 
matter that have been little explored, including understanding the physical interactions of colloids 
and microswimmers with their viscous 3D environment and the interplay between these stimuli and 
the molecular response of biological systems. 


INTRODUCTION 


Probing biology and physics at high resolution and speed in 3D is a tremendous experimental challenge that has 
inspired significant advances in microscopy techniques across time and length scales. Fluorescence imaging is one 
of the most widely used techniques, with many existing variants addressing different experimental needs from the 
nano- to macro-scale. Fundamentally, all fluorescence methods are subject to the “triangle of frustration”, where 
photobleaching, phototoxicity and fluorophore photophysics limit the speed, duration, and total light dose [1]. 

Recent advances in quantitative phase imaging (QPI) techniques have positioned these as an alternative to fluores- 
cence in many domains [2]. In contrast to fluorescence microscopy, image contrast in QPI originates from the sample’s 
refractive index (RI), avoiding the aforementioned low-signal levels, photobleaching, and the triangle of frustration 
associated with fluorescent labeling. When morphology and dynamics are of primary interest, QPI approaches are 
becoming a method of choice for long-term, volumetric imaging at few-hertz rates [3]. 

For acquisition rates beyond a few hertz, a variety of single-shot QPI approaches enable fast “volumetric” imaging, 
including digital in-line holography [4], off-axis holography [5], and iSCAT [6] (table I). However, single-shot QPI 
methods can only obtain unambiguous volumetric reconstructions by putting strong priors on the geometry of particles 
that are imaged, e.g., spheres, rods, or other particles with a known scattering model. Such approaches commonly 
fail in dense samples and have limited axial resolution [7]. 

In contrast, optical diffraction tomography (ODT) [8-12], infers the 3D RI of a sample by combining multiple 
views obtained by projecting coherent light through the sample at different angles. At a high level, ODT is closely 
related to synthetic aperture microscopy [13] and Fourier ptychography [14]. ODT distinguishes itself from other 
“volumetric” QPI approachs in its flexibility to probe arbitrary refractive indices and non-sparse samples and obtain 
true volumetric images. Therefore, ODT stands out as the method of choice for imaging dense samples in complex 
environments, as evidenced by a number of recent biological applications including measuring cell dry mass [15], flow 
cytometry [16], traction force microscopy [17], and 3D histopathology [18]. 

Due to the need to acquire 2 100 views to perform high-quality tomographic reconstruction, ODT imaging is 
typically two orders of magnitude slower than single-shot techniques and achieving high-quality RI reconstructions 
at kilohertz volumetric imaging speeds is an open challenge. The primary hurdle to improve ODT acquisition speed 
is the rate of pattern generation for individual ODT views. Commonly, the angle that the illumination traverses the 
sample is controlled by galvanometric mirrors [19] or spatial light modulators (SLM’s) [20-24]. The fastest and most 
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common SLM choice is the digital micromirror device (DMD), a binary device often run in a time-multiplexed gray- 
scale mode that limits the pattern refresh rate to 1 kHz [25, 26]. DMD’s are capable of pattern changes at >10 kHz 
in binary mode, but this introduces stray diffraction orders that must be filtered out using static [21] or dynamic 
(16, 27] masks. Even working at the full DMD display rate, achieving high-quality ODT using >100 patterns limits 
the volumetric frame rate to <100 Hz, still short of the volumetric kilohertz target. 


Method speed |limit contrast |true 3D |ref 
digital holographic microscopy |5 MHz |camera RI no [28] 
iSCAT 20 kHz |camera RI no [29] 

SLIM 1kHz |camera/multiplexing F yes [30] 
eventLFM 1kHz |sample/accumulation time|F yes [31] 

ODT 0.1kHz| DMD/camera RI yes see text 
FS-ODT 1kHz |DMD/multiplexing RI yes this work 


TABLE I. Comparison of kHz-scale volumetric imaging methods. RI and F stand for refractive index and fluorescence 
respectively. The quantitative phase imaging techniques are discussed in the text. SLIM and eventLFM are light field microscopy 
techniques. 


To address the challenge of realizing high-resolution, true-3D, kilohertz-rate, and label-free volumetric imaging, we 
developed an alternative strategy to rationally designing and projecting ODT illumination patterns, Fourier synthesis 
optical diffraction tomography (FS-ODT). FS-ODT builds on our previous DMD pattern generation tools [32] to 
enable multiplexing many ODT views into a single image. We refer to this strategy as Fourier synthesis because 
we construct the complex illuminations patterns required for multiplexing by combining simple atomic patterns in 
the Fourier plane. Because FS-ODT places the DMD in a conjugate Fourier plane to the objective and achieves 
position and angle control over the patterns using a spatial carrier frequency, it is possible to increase the information 
content in a single image orders of magnitude further than multiplexing Lee holograms [16]. While multiplexed 
intensity diffraction tomography (IDT) also synthesizes complex patterns using multiple LED point sources in the 
far-field [33, 34], the use of a fast pattern modulator enables additional control over the beam phase and position and 
allows us to achieve two orders of magnitude faster volume acquisition. 

Increasing the information content per image by multiplexing introduce a more ill-posed reconstruction problem. 
To address the ill-posed problem of multiplexed patterns, we designed an iterative reconstruction approach based on 
multi-slice beam propagation and accelerated proximal gradient descent. We additionally created a “demultiplexing” 
approach to generate a high-quality initial 3D RI distribution and created new approaches to avoid non-optimal 
solutions due to phase wrapping. 

We demonstrate the capabilities of FS-ODT for a variety of samples and scenarios, illustrating its unique combina- 
tion of high-quality RI reconstruction and high-speed volumetric imaging capability. We first profile the reconstruction 
quality versus the degree of angle multiplexing by imaging a sample of known composition, a 101m polystyrene mi- 
crosphere (PMS) and a biological sample of known morphology, Tetrahymena. Next, we consider dynamic samples, 
including diffusing PMS and bacteria, and demonstrate FS-ODT tracking of colloidal particles and extraction of their 
hydrodynamic properties. Finally, we combine position and angle multiplexing at the fastest FS-ODT imaging rate 
possible with our hardware to measure diffusing microspheres at kilohertz volumetric frame rates. Taken together, 
the results presented here demonstrate that FS-ODT is a powerful new approach for the dynamics of cells and colloids 
in complex environments. 


RESULTS 
Fourier synthesis of ODT patterns 


In FS-ODT, we generate patterns using a DMD conjugate to the objective back focal (or Fourier) plane, as shown 
in Fig. la. Plane waves are encoded by circular “spot” patterns on the DMD. Varying the beam position in the sample 
corresponds to a adding a phase ramp in the Fourier plane, and so to make the beam position tuneable we align the 
system so that beams which diffract from a given spatial carrier frequency on the DMD are centered in the optical 
system. Beam angle through the sample is encoded by spot position, so multiplexing tens to hundreds of beams at 
different angles is possible by displaying multiple spatially separated spot patterns on the DMD. Further, by working 
in the Fourier plane we avoid the need to filter stray diffraction orders as these are mitigated by Fourier broadening 
due to the small spot sizes. We have carefully modelled this approach by extending the DMD simulation tools we 
developed previously for multicolor structured illumination microscopy [32] (section $2). 


FIG. 1. Fourier synthesis optical diffraction tomography. a. In FS-ODT, a spot-like pattern with a superimposed 
spatial carrier frequency is displayed on the DMD (lower left), which is conjugate to the back-focal plane of the objectives 
(purple line). The light which diffracts from the — mirrors (red) at the carrier frequency is transmitted, while light which 
diffracts from the + mirrors (blue) is blocked. The position and carrier frequency of the spot pattern control the angle and 
position of the ODT illumination respectively. b. Fixed Tetrahymena imaged with FS-ODT shown in a single axial plane. 
Internal cell structures are resolved (left), and ~200 nm cilia are at the edge of our detection ability (right). c. MIP’s of the 
Tetrahymena RI in the xy- (upper left), xz- (lower left), and zy-planes (right). Internal structure of the cell can be resolved, 
including the nucleus and the oral apparatus. The unusual morphology suggests this cell may be in the process of dividing. 
Scale bar 20 pm. 


To generate a single plane wave at spatial frequency f = (fz, fy) displaced from the center of the field of view by 
dr, we generate a spot pattern at location r, and spatial frequency f, on the DMD. The object space parameters are 
determined by 


fir 

br = 57 (fp ~ fe) (1) 
M 

f= Fa (fp — Ye) (2) 


where f, is the spatial carrier frequency, r, is the position on the DMD aligned with the optical axis, f; is the focal 
length of the objective lens, M is the magnification between the DMD and the lens back focal plane, and A is the 
wavelength of light. This approach is power inefficient as the incident beam illuminates the full DMD chip, but 
spot patterns of typical radius 10 mirrors diffract only a small portion of the light. Nevertheless, since ODT detects 
transmitted light, there is sufficient signal to image with <100 ps exposure times (section S2 A). 

To validate that the proposed FS-ODT optical design is capable of high-quality 3D RI reconstruction, we first 
performed non-multiplexed ODT on a complex 3D sample, fixed Tetrahymena cells (Fig. 1b,c). After RI reconstruc- 
tion, we resolve the internal structure of the cell, including the nucleus, the oral apparatus, and the cilia. The cilia 
are expected to be 0.2m in diameter, which is significantly smaller than the Abbe limit for our detection optics, 
2NA/A ~ 0.39pm. As such, the measured refractive index contrast is small due to averaging of the structure and 
background over the resolution voxel size. Reconstruction parameters are provided in table II. 

We infer the RI distribution using custom GPU-accelerated Python code which implements accelerated proximal 
gradient descent using the fast iterative shrinkage-thresholding algorithm (FISTA) [35] (sections A and $3). We in- 
clude the physics of light propagation through the RI using either the beam-propagation model (BPM) or the split-step 
non-paraxial (SSNP) forward models [34, 36]. FISTA and related proximal approaches are particularly powerful for 


FIG. 2. Multiplexed illumination with FS-ODT. a. Exemplary position multiplexing pattern. Five spot patterns are 
displayed across the DMD face (top left) with different carrier frequencies (inset, bottom left). When relayed to the objective 
back focal plane (right), the resulting beams have different incidence angles, and the lens transforms these angles to different 
positions in the focal plane. Only three beams are shown for clarity. b. Exemplary angle multiplexing pattern. 12 spot patterns 
(top left) are displayed at different spatial positions but with the same carrier frequency (inset, bottom left). When relayed 
to the objective back focal plane (right), the resulting beams have different incidence positions, and the lens transforms these 
positions to different angles in the focal plane. Only five beams are shown for clarity. 


solving ill-posed inverse image reconstruction problems because they provide a framework for applying regularization. 
Regularization stabilizes the reconstruction and favors physical solutions from all possible RI distributions consistent 
with the data. We use total variation regularization (TV), which promotes smoothness, and ¢; regularization, which 
promotes sparsity. Additionally, we impose physical constraints, typically that the imaginary part of the RI is strictly 
absorptive and the real part of the RI is greater than the background index. 

Next, we apply the FS-ODT pattern generation strategy to generate multiplexed illumination patterns. Since 
FS-ODT uses spatially separated spot patterns to generate plane waves at different angles, we generate multiplexed 
patterns by displaying spot patterns at different positions on the DMD. We synthesize composite patterns from a 
set of beam frequency and position shift pairs {(f1,6r1),...(f{v,6rn)}. Each pattern generates N beams passing 
through the sample plane in different positions at different angles (i.e. spatial frequencies). Unlike in multiplexed 
IDT approaches, these beams are mutually coherent, and so there combination generates a complex illumination 
pattern at the sample. Our approach enables two forms of multiplexing: (i) position multiplexing to extend the 
system field of view and (ii) angle multiplexing, which illuminates a single sample region with overlapping beams 
at different angles. We demonstrate these two strategies in Fig. 2. A unique benefit of FS-ODT is that these two 
strategies can be combined to enable angle multiplexing over a large field of view. 


Low-resolution Rytov demultiplexing 


Iterative ODT reconstruction, as described above, relies on a high-quality initial guess for the RI distribution. 
Since the problem is not globally convex, the iteration scheme may not converge if the guess is too far from the true 
RI distribution. The initial guess is most commonly generated from linear weak scattering approximations, such as 


FIG. 3. Low-resolution Rytov demultiplexing. a. In the Fourier transform of each off-axis hologram, we select regions 
around each plane wave frequency where the scattered light is dominated by that individual plane wave (yellow circle). The 
radius of these regions is significantly smaller than the bandpass frequency of the microscope, fp, so these regions contain 
low-resolution information about the RI object. b. The information from each low-resolution region is mapped to the correct 
position (yellow arcs) in the 3D Fourier space representation of the scattering potential. These regions are centered about zero 
frequency (black dot). For reference we show arcs corresponding to a straight-on incident beam (purple) and beams at the band 
pass frequency (green). c. We obtain an initial guess for the RI after inverse Fourier transforming the scattering potential. For 
the 10 ym diameter PMS shown here, the low-resolution demultiplexed guess in the Rytov approximation provides an excellent 
starting point for further optimization. 


the Born and Rytov approximations, that can be efficiently computed from the electric field data. For plane wave 
illumination, these approximations relate the 2D Fourier transform of the image electric field to a spherical cap in the 
3D Fourier transform of the scattering potential (section $3 A). 

Multiplexing many patterns results in a more ill-posed inverse problem than single-beam ODT due to the additional 
need to unmix scattering contributions from different incident plane waves. When the illumination pattern contains 
multiple plane waves, it is not obvious which spherical cap each image Fourier frequency should be assigned. This 
ambiguity can be removed experimentally using phase shifting, but this requires additional images and negates the 
speed advantage of multiplexing. 

To improve the convergence of our reconstruction algorithm, we develop an approach to initialize the RI with a 
high-quality guess in the presence of multiplexing, which we refer to as low-resolution Rytov demultiplexing. In 
this scheme, we observe that in many common experimental situations, information from the scattering of ith plane 
wave dominates the hologram in the region around the plane wave carrier frequency. We expect this situation to 
prevail if scattering is weak enough and, for example, the sample RI is greater than the background and its power 
spectrum decays with increasing spatial frequency. For each plane wave, we select the region in Fourier space which 
is closer to that carrier frequency than to any other. We take the resulting N lower-resolution demultiplexed fields 
and compute the corresponding Rytov phases (eq. SS18). Finally, we generate an initial RI distribution by using the 
linear scattering model to map this information to the scattering potential at the correct position in Fourier space. 
An example of this process is shown in Fig. 3a, b for 10x multiplexing for a sample consisting of 101m diameter 
PMS. The resulting Rytov demultiplexed RI (Fig. 3c) captures the features of the sample and is a good initial guess 
for the solver. 

While we find Rytov demultiplexing produces a better initial guess than e.g. beginning with a constant RI, it 
eventually fails because as the degree of multiplexing increases, the resolution of the demultiplexed images decreases. 
For small objects and high degrees of multiplexing, the demultiplexed images do not contain sufficient detail to 
accurately represent the RI. This problem is particularly acute for thick objects that produce phase-wrapped electric 
fields, because the phase wraps must be detected by comparing phase changes across the object. When rapid phase 
variations are not captured by the low-resolution demultiplexed images, phase unwrapping breaks down and phase 
errors in the initial guess can lead to rapid deterioration of the RI reconstruction. 

In some cases, the sensitivity of the optimization problem to correct phase-unwrapping can be mitigated by using 
a loss function which includes a phase-sensitive and a phase-insensitive term (eq. $5). This illustrates that, while the 
phase information obtained by ODT produces a less ill-poised reconstruction problem it simultaneously produces a 
more complex landscape for the loss function which poses significant challenges for convex optimization approaches. 


FS-ODT validation 
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FIG. 4. 10um microsphere refractive index reconstruction versus angle multiplexing. a. RI reconstruction using 
the SSNP and an initial guess from the low-resolution demultiplexed Rytov approach described in the text. MIP’s through 
the microsphere are shown in the xy-plane (upper left), yz-plane (upper right), and xz-plane (lower left). Results for angle 
multiplexing factors of 1, 3, 10, and 19 are shown (left to right). Scale bar 107m b. Exemplary raw hologram images which 
are used to recover the bead RI. Scale bar 10m. For higher degrees of multiplexing, these become more complex due to the 
mutual interference of many plane waves. c. Exemplary DMD patterns used to generate the multiplexed illumination patterns 
in B. Only the portion of the DMD within the detection objective pupil radius (grey circle) is shown. The DMD is slightly 
smaller (~84%) than the pupil along the vertical direction. 


As an initial validation of our approach, we collect FS-ODT images of a sample of known structure and RI, a 
10m diameter PMS suspended in immersion oil, using a range of different multiplexing conditions. We consider 
multiplexing by factor of 1, 3, 10, 19 using a fixed set of N = 147 plane wave directions. Where the total number of 
angles is not divisible by the multiplexing factor, we include some plane waves twice to simplify the reconstruction. 

While FS-ODT can generate arbitrary multiplexed patterns, reconstruction quality is improved using multiplexed 
patterns that are easier to unmix by ensuring the beam frequencies in each image are as far apart as possible. This 
maximizes the resolution achieved in the Rytov demultiplexing initialization. We design our patterns using an iterative 
algorithm to ensure this is the case. We initialize each M-fold multiplexed pattern with M beam angles by iteratively 
selecting the beam angle that maximizes the distance from those already chosen. Next, we iterate over all choices of 
two patterns and randomly swap beam angles. We define a loss function, the average distance between beam angles 
on the DMD up to a certain maximum value. If swapping the angles increases the loss function, we keep the swap. 
We obtain high-quality multiplexed patterns after performing 5 iterations of 300 swaps. 


Due to the thickness and large RI contrast between the PMS and the immersion oil, npyyg ~ 1.57 and no ~ 1.515, we 
use the SSNP for our reconstruction [36] to both allow for multiple scattering and accurately model the propagation 
of the complex illumination beam. The BPM with an obliquity factor correction achieves similar performance with 
reduced memory and computational requirements for single plane wave illumination [34], but in our case the multiple 
incident plane waves do not share a common obliquity factor, so the SSNP is more appropriate. 

We find that FS-ODT achieves high-quality 3D RI reconstructions of PMS for 1x (no multiplexing), 3x, 10x, 
and 19x multiplexing (fig. 4). With no multiplexing, we recover a nearly spherical object of the correct diameter 
with n = 1.57, as expected because this problem is the least ill-posed and has the largest effective SNR per beam. 
Since all raw images are acquired with similar peak intensity, the effective SNR per beam is reduced as the degree 
of multiplexing increases. As expected, due to the missing cone problem, the sphere appears somewhat elongated in 
the axial direction. The 3x, 10x, and 19x multiplexed data also produces high-quality reconstructions, but the RI is 
not distributed as uniformly along the optical axis. This reconstruction artifact is reminiscent of missing-cone effects, 
and could be a result of effective loss of SNR for individual plane waves in the highly-multiplexed beam. 

We numerically explore and further validate our multiplexed FS-ODT approach by reconstructing synthetic images 
generated using custom GPU accelerated Mie scattering code (section $1). We find that multiplexed patterns signif- 
icantly enhance the RI reconstruction compared with non-multiplexed patterns when the number of raw tomograms 
is kept fixed. 


Hindered diffusion and microswimmer motility 
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FIG. 5. Hindered diffusion of colloidal microspheres. a. Colloidal particles interact with nearby boundaries when the 
boundary influences the fluid flow fields generated by the particles’ motion. b. 1 jm diameter PMS diffusing in water-glycerol 
mixture imaged with FS-ODT. MIP’s of the RI are shown in three orthogonal planes. Circles indicate bead location determined 
by tracking algorithm. Color represents axial position from near the coverslip (red) to up in the sample (purple). Scale bar 
20um. c. Hydrodynamic interactions between microspheres and the coverslip induce a spatially varying diffusion coefficient 
(top). The microspheres are buoyant in this solution, leading to a spatially varying density profile (bottom). d. FS-ODT 
resolves individual microspheres in a solution with x10 greater PMS concentration compared with b. Scale bar 20 pm. 


Having established FS-ODT produces high quality RI reconstructions in static samples, we now apply our approach 
to study hindered diffusion, hydrodynamic interactions, and microswimmer motility. Phase imaging approaches offer 
significant advantages over fluorescence imaging in lower phototoxicity, faster frame rates, and axial position sensitivity. 
As such, phase methods are ideally suited to studying motion in 3D environments, such as the motility of colloids 
or cells in complex environments. One poorly understood question is how cells swim in viscoelastic environments 
created by the presence of inert colloids or polymers in solution [37]. This topic is of broad interest because these 
fluids mimic the natural environment of bacterial species better than typical in vitro experiments and may provide 
insight into behavior and evolutionary pathways which is currently lacking. Complex environments may also impact 
collective bacterial dynamics, including swarming and biofilm formation, which are interesting from a fundamental 
perspective of better understanding active matter as well as from a public health perspective [38]. 3D QPI approaches 
offer significant advantages over 2D-only or 2D with limited z-range tracking, which only measure particles at or near 
the microscope’s focal plane, particularly in dense environments like bacterial suspensions in complex fluids. As a first 
step towards addressing these questions, we demonstrate that we can image and distinguish diffusing microspheres 
and swimming FE. coli cells in 3D using FS-ODT. 

Imaging dense samples introduces new reconstruction challenges, as our algorithms require knowledge of the exci- 
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FIG. 6. E. coli motility. a. RI reconstruction of bacteria and tracer particles MIP along three orthogonal axes, XY (upper 
left), YZ (right), and XZ (bottom). The bacteria and tracers can be distinguished based on their RI and morphology. Two 
bacteria (upper right) and ~20 tracer particles are visible in the XY projection. b. Bacteria and tracer particles 30 ms later. 
One bacteria has rotated, and the tracers have changed position due to diffusion and advection. c. During bacterial swimming 
we define the body axis (blue) and the average velocity (red). The wobble angle, 0, is the average angle between the two. d. 
Speed (top) and 2D orientation of the bacteria axis (red and gold) and the bacteria velocity (blue and purple). Different traces 
represent the two swimming bacteria observed in the experiment. e. Histogram of measured wobble angles. Colors correspond 
with bacteria identified in d. 


tation electric field with no sample present to distinguish electric field patterns generated by the complex illumination 
pattern from those due to interaction with RI objects. One standard solution is to acquire a background image taken 
at a spatial position where no RI objects are present. However, for dense diffusing objects, there are no sample-free 
regions. Instead, we rely on the time-average image, assuming that over the long term, the RI inhomogeneities average 
out. For relatively sparse samples, this is a good approximation. For denser samples, it may be necessary to account 
for an average background RI different from the fluid. Other approaches, including joint inference of the incident 
electric field and the sample RI, are possible [39]. 

We applied FS-ODT to study hindered diffusion of 1 ym beads in a water glycerol mixture and quantified the bead 
dynamics by determining the diffusion coefficient as a function of distance to the coverslip using 3D tracking (Fig. 5, 
Videos 1 and 2). The diffusion coefficient is sensitive to hydrodynamic forces on the beads and serves as an indicator 
of hydrodynamic interactions. Colloidal particles experience hydrodynamic interactions with boundaries because they 
generate fluid flow fields as they move, which must satisfy the Stokes equation with appropriate boundary conditions. 
When a sphere is within a few radii, R, of a boundary, the no-slip boundary condition at the wall modifies the 
flow fields (Fig. 5a) which affects the mobility of the sphere and, therefore, the translational and rotational diffusion 
coefficients. This effect, known as hindered diffusion, is perturbative in the ratio of the sphere’s radius to the distance 


from the wall R/h. To leading order, the diffusion coefficient parallel to the wall is [40] 


9R R? 
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where D, is the diffusion coefficient in a unbounded fluid. Our measured diffusion coefficients versus height above the 
coverslip (Fig. 5b) match well this functional form. There are small deviations from the expected curvature, which 
we attribute to the fact that the microspheres sample different heights during our measurement, leading to some 
broadening of the measured curve. 

Hindered diffusion has previously been studied using evanescent wave dynamical light scattering [41] and optical 
tweezer experiments combined with in-line holography [42]. FS-ODT provides an alternative approach which can 
resolve hydrodynamic forces over tens of 1m axially without the need to perturb the sample, scan tweezer positions, 
or combine multiple imaging techniques. ODT’s unique ability to image across a 3D field of view and in multiple 
scattering environments enables simultaneously tracking of many particles in dense colloidal suspensions, which is 
extremely challenging using these other approaches. 

Next, we studied EF. coli motility by imaging bacteria in solutions seeded with 0.5m PMS tracer particles at 
volumetric rates of 143 Hz over ~7s (Fig. 6, Videos 3 and 4). We found that the bacteria and tracer particles can be 
distinguished based on RI and morphology. We observe diffusion of the tracer particles and directed swimming of 2 
E. coli cells in this dataset. To quantify the swimming behavior of the EF. coli, we trained a classifier to segment the 
cell bodies and then tracked their position and orientation (Fig. 6D). From these 4D data, we computed the wobble 
angle, which is the average angle between the cell’s average velocity and body axis (Fig. 6c). Wobble angles have 
previously been used to characterize E. coli swimming in complex fluids [37]. This experiment demonstrates FS-ODT 
can effectively identify bacteria swimming in 3D and be used to compute wobble angles, which we anticipate will 
enable larger-scale studies. 


Diffusing microspheres imaged at kilohertz volumetric rates 


We next demonstrate that FS-ODT can reliably image diffusing colloidal particles at kilohertz volumetric rates. 
Specifically, we image 1m PMS diffusing in water over a period of 0.969s at a volumetric rate of 1.032kHz with a 
field-of-view of 20.5 x 96 x 104 pm (Videos 5 and 6). At each time point we collect 8 images using 19x angle multiplexed 
patterns. Additionally, we expand the field-of-view of the imaging system by expanding the illumination pattern 2 in 
each direction by including spots with 4 different carrier frequencies. This combined angle and position multiplexing 
uses a total of 608 plane wave patterns. We show an exemplary angle- and position-multiplexed illumination pattern 
in Fig. 7a-c. 

We performed RI reconstruction at each time point (Fig. 7f) and tracked the diffusing PMS to extract their 
diffusion coefficients. As for other time-resolved experiments, we construct a reference electric field by time-averaging 
the electric fields obtained from off-axis holography. We quantify the diffusive motion by calculating cumulative 
distribution functions (CDF’s) of the step-size distribution along each axis (Fig. 7d) and the mean-square displacement 
(MSD) of the beads versus lag-time (Fig. 7e). For short lag-times, the CDF’s can be strongly affected by localization 
error. Therefore we examine the CDF’s for a lag-time of 0.194s (200 frames) and find that the CDF’s closely match 
the expected Gaussian form. We extract the ensemble-average diffusion coefficient by performing a linear fit to the 
MSD and determine D = 0.68 :m?/s. 

We display volumetric data in Fig. 7 using maximum intensity projections (MIP’s) along three orthogonal axes. We 
can easily resolve the diffusing PMS, and we find that they are distributed over an axial range of about 10m, which 
is the approximate height of the sample chamber (Fig. 7f). We also imaged PMS diffusing in a 120 pm thick chamber, 
focusing on the region within 40 pm of the coverslip (Fig. 7g). These images provide an ideal starting point for particle 
tracking, particle imaging velocimetry, or other approaches to either explore the physics of colloidal suspensions or 
use these colloids as tracer particles to report on more complex biological systems. 


DISCUSSION 


Optical diffraction tomography is a powerful label-free 3D imaging approach with broad applications in biomedical 
imaging, biophysics, and soft matter physics. However, most efforts focus on improving ODT for applications specific 
to biomedical imaging. Biophysics and soft matter systems have very different imaging requirements, motivating our 
development of multiplexed FS-ODT pattern generation and reconstruction approaches. 

We demonstrated that FS-ODT provides high-quality RI reconstructions of a variety of samples, including cells, 
10m microspheres, diffusing colloids, and swimming bacteria. We present multiplexed illumination strategies to 
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FIG. 7. Diffusing microspheres at kilohertz volumetric rates with an extended field-of-view. a. Intensity 
profile of the beam after interacting with microspheres using 19x angle multiplexing and 4x position multiplexing (circles). 
Position multiplexing expands the usable field of view by ~ 2x b. Fourier transform of the intensity profile in a showing 76 
distinct beams in Fourier space. c. DMD pattern used to generate the illumination shown in a and b d. Step-size cumulative 
distribution function (points), Gaussian model fits (lines), and residuals (bottom) for steps separated by 200 frames along the 
x (red-orange), y (yellow), and z (purple) directions. e. Ensemble-average mean-square displacement (MSD) versus lag-time 
(grey) and a linear fit (red) up to a lag-time of 0.144s. f. RI reconstruction at a single time point for a ~10 pm tall chamber. 
RI values are shown between 1.352 and 1.413. PMS are localized (circles) and tracked. g. RI reconstruction at a single time 
point for a denser sample in a 120 pm tall chamber. RI scale the same as f. 


increase the number of illumination angles in a single image and expand the system’s field of view. Combined, both 
multiplexing strategies enable kHz-scale ODT over extended lateral fields view. 
While the kHz-scale volumetric imaging is fast enough for many applications, FS-ODT imaging could be further 
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accelerated by adopting alternative pattern generation hardware. For example, using faster DMD models that have 
pattern update rates up to 30kHz could immediately accelerate the volumetric imaging rate by a factor of 3. Al- 
ternatively, acousto-optic deflectors or electro-optic modulators could accelerate pattern generation by 1—2 orders 
of magnitude, at the cost of complicating the off-axis holography reference generation. These fast speeds might be 
necessary for studying swimming in fast ciliates, which can move at speeds of up to 4000 m/s [43] or studying the 
diffusion of single proteins [29]. 

While FS-ODT already produces high-quality reconstructions in most circumstances, further improvements are 
possible by applying FISTA or a deep learning approach to unmix the off-axis holograms, denoise raw data, suppress 
coherent speckle, or perform phase unwrapping. Supervised learning approaches may improve the speed or quality of 
multiplexed reconstructions [16, 44]. 

Alternatively, emerging deep-learning regularization approaches could be incorporated to improve FS-ODT recon- 
structions further. For example, DNN denoisers can be incorporated as priors in our existing FISTA framework using 
the plug-and-play prior or regularization by denoising approaches [45]. Alternatively, recent self-supervised approaches 
have achieved impressive results while preserving a physics-based forward model by using DNN’s to parameterize the 
RI distribution, either as a pixelated image as in deep image prior approaches [46] or as a function of coordinates, as 
in neural field based approaches [47]. To account for motion blur, it may be possible to adapt neural field approaches 
to create a space-time model and infer dynamics during hologram acquisition [48]. 

We anticipate that FS-ODT will unlock challenging imaging regimes that have been difficult or impossible to explore. 
These include probing the ~100 Hz mechanical motion and hydrodynamics of cilia during swimming of protists such 
as Tetrahymena or exploring the behavior of microswimmers and active colloids in dense or complex 3D environments. 
It also extends the frontier of QPI to complex object tracking in dense environments and the exploration of highly 
dynamic systems. 
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sample figure M|Np| NM rate exp |model}n, f | tw | 7, |max n”|batch dz x dxy FOV (zx y x 2) Te coe 
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10pm PMS 1x |Fig. 4 147/147 1 400 ps|SSNP 11.515) 1 | 0.1 | 0.003 0 4 0.1 x 0.1 pm 25.1 x 70x 70pm = |20 min 
10pm PMS 8x |Fig. 4 49 |147 1 300 ps|SSNP 11.515) 1 | 0.1 | 0.003 0 4 0.1 x 0.1 pm 25.1 x 70 x 70pm {16min 
10 ym PMS 10x |Fig. 4 15 |150 1 150 ps|SSNP ]1.515] 0.1 | 0.1 0 0 2 0.1 x 0.1 pm 25.1 x 70x 70pm {32 min 
10pm PMS 19x |Fig. 4 8 |152 1 75s |SSNP /1.515/0.01] 0.1 0 0 2 0.1 x 0.1 pm 25.1 x 70x 70pm {32min 
aa 
1pm PMS ees 11] 11 ]10000} 6.05Hz | 3ms |Rytov|1.407) 1 | 0 | 0 coo | 1 | 1.75x0.376pm| 54x 67x 79pm | 14h 
lym PMS oo 11] 11 | 2000 | 6.05Hz | 3ms |Rytov|1.407) 1 | 0 | 0 coo | 1 | 1.75x0376pm] ~54 x 67x 79m 
Fig. 6c,d 
0.5m PMS . p : : ‘ 
; Video 3 | 11) 11 | 1000 | 143Hz |600ps]BPM |1.333) 1 |0.03] 0.01 oo 11 {0.435 x 0.196 pm} 25.6 x 25.3 x 26.5 pm|49 min 
& E. coli . 
Video 4 
lpm PMS oe 8 |608] 1000 |1.032kHz} 75 ps |BPM /1.333} 1 |0.01]/0.0003 0 3 0.5 x 0.1 pm 40.5 x 96 x 104m | 48h 
1pm PMS ea 8 |608| 1000 |1.032kHz| 75 ps |BPM |1.333) 1 | 1 | 0.01 0 3 0.5x 01pm | 70.5x96x104pm | 66h 
10pm Mie 10x |Fig. Sla | 15 |150 1 SSNP ]1.515) 1 1 0 0 1 0.1 x 0.1 pm 25.1 x 70 x 65 ym 
10pm Mie 1x Fig. Slb | 15] 15 1 SSNP ]1.515) 1 1 0 0 1 0.1 x 0.1 pm 25.1 x 70 x 65pm 


TABLE II. ODT reconstruction parameters. Here we provide details about the ODT patterns and reconstruction param- 
eters used in all data sets presented above. The pattern parameters are M, the number of images, N,, the number of plane 
wave patterns, and N;, the number of time lapse points collected. Rate is the volumetric image acquisition rate and exp is the 
exposure time for each off-axis hologram. The reconstruction parameters are the fraction of electric field loss function, f, the 
regularization parameters Ty and 7¢,, the maximum allowed value for the imaginary part of the RI, max n”, the voxel size of 
the reconstruction grid, dz x dxy, and the reconstruction field of view, FOV. Trecon is the reconstruction time. 


MATERIALS AND METHODS 
A. Computer control 


The microscope is controlled using a custom computer package [51], and a GUI based on Napari-MicroManager, 
which relies on the MicroManager device drivers but replaces the Java GUI with Napari, https://github.com/ 
QI2lab/napari-micromanager. The control computer is a Lenovo ThinkStation P620 running Microsoft Windows 
10 Enterprise. This computer has an Ryzen Threadripper Pro 3945WX CPU (AMD) with 12 cores, 128 GB RAM, and 
an GeForce RTX 3090 GPU (NVidia) with 24GB of memory. ODT patterns are pre-loaded on the DMD firmware, and 
the entire microscope acquisition is hardware triggered by a National Instruments DAQ (PCIe-6343). The Phantom 
camera includes 72 GB on-board memory and the 10 Gigabit ethernet option. It communicates with the PC using an 
X540-T2 10GbE card. A full 72GB camera acquisition can be transferred to the computer in ~2 min. 


B. Refractive index reconstruction 


FS-ODT reconstructions were carried out using custom Python code [51] running on Python 3.11.4 or 3.11.5 and 
CUDA toolkit 11.8.0. We implemented light scattering forward models and proximal gradient algorithms, harness- 
ing CuPy for GPU acceleration and Dask for parallelization and orchestration. We rely on the RAPIDS cuCIM 
implementation of the total variation proximal operator, which is modeled on the scikit-image implementation, 
denoise_tv_chambolle. For weighted phase unwrapping during computing the Rytov field, we adapted the approach 
of [52] to run on the GPU. 

Reconstructions were performed on either the experimental control computer described above, or a custom-built 
computer running Ubuntu 20.04.6 LTS with an X99S SLI Plus motherboard (MSI), an i7-5820K CPU (Intel), 125 GB 
RAM, and an RTX A6000 GPU (NVidia) with 48 GB of memory. 


C. Preparation of Tetrahymena samples 


Tetrahymena were cultured at room temperature in media containing 2% proteose peptone, 0.2% yeast extract, 
0.012mM FeCls, 0.2% glucose, 100 unit/mL penicillin, 100 mg/mL streptomycin, and 0.25mg/mL Amphotericin B. 
Tetrahymena were passaged to fresh media every 3-4 days. Tetrahymena were stained as described previously at room 
temperature unless noted [53]. Briefly, ~2 x 10° mid-log cells were washed with 10mM Tris pH 7.5, permeabilized 
in 0.25% TX-100 in PHEM for 30s (60mM PIPES, 25mM HEPES, 10mM EGTA, 4mM MgSO,), fixed in 1% 
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paraformaldehyde in PHEM for 15min, blocked in 3% BSA in PBS-T for 30min (0.01% Tween-20, 130mM NaCl, 
2mM KCl, 8mM NagHPO,, 2mM KH2PO,, 10mM EGTA, 2mM MgCl, pH 7.5), stained with primary antibodies 
diluted in 3% BSA in PBS-T at 4°C overnight (glutamylated tubulin 1:1000, GT335, Adipogen, AG-20B-0020-C100 
and centrin 1:1000, 20H5, MilliporeSigma, 04-1624), stained with secondary antibodies diluted in 3% BSA in PBS-T 
for 1 hour (1:1000 goat anti-mouse IgG2a Alexa Fluor 488, Invitrogen, A-21131 and 1:1000 goat anti-mouse IgG1 
Rhoadmine Red X, Jackson Immuno, 115-005-205), and counterstained with 1mg/mL DAPI for 5min . Samples were 
washed with PBS three times after fixation, primary antibody incubation, and secondary antibody incubation. Each 
wash lasted 5 min and all centrifugations were carried out for 5min at 250 x g in a swinging bucket centrifuge. 

For the samples in Fig. 1b, we placed ~100 pL of fixed Tetrahymena cells in PBS in a 40mm diameter cell culture 
dish (FluoroDish FD5040). For the samples in Fig. 1c, 10 pL of fixed Tetrahymena cells in PBS were placed on a round 
#1.5 coverslip of diameter 40mm. A square #1 coverslip with 25mm sides was dropped on top and the chamber was 
sealed with epoxy (Devcon 2 Ton Epoxy, GLU-735.90). The electric field data was binned by a factor of 2 to reduce 
the memory required during reconstruction. 


D. 10pm polystyrene microspheres 


PMS of diameter 10m (Thermofish Fluospheres F8836) were sonicated for 20 min and then diluted by a factor of 
10 in immersion oil (1.04699.0500 MilliporeSigma) with RI in the range 1.515-1.517. We spread 20 pL of the resulting 
emulsion on a #1.5 coverslip of diameter 40mm with a pipette tip and left it uncovered overnight for the water to 
evaporate. Then we dropped a #1 square 25 x 25mm coverslip on top and sealed the chamber with epoxy (Devcon 2 
Ton Epoxy, GLU-735.90). After the epoxy dried, the sample was placed on the microscope, and a few drops of water 
were placed on the top coverslip to facilitate imaging with a water dipping objective. 


E. Preparation of diffusing microspheres in a water-glycerol mixture 


PMS of diameter 2R = 1 pm (ThermoFisher Fluospheres F8823) of weight/volume 0.02 gmL~' were first sonicated 
for 20min and then diluted by a factor of 10 in MQ water. 10 pL of this dilution was combined with 40 pL of MQ 
water and 50 pL of glycerol. 10 pL of this mixture was placed on a round #1.5 coverslip of diameter 40mm. A square 
#1 coverslip with 25mm sides was dropped on top, and the chamber was sealed with epoxy (Devcon 2 Ton Epoxy, 
GLU-735.90). 

This experiment was performed on an earlier version of the apparatus. The detection objective was a 50x long- 
working distance air objective with NA = 0.55 (Mitutoyu 378-805-3) and 200mm focal length tube lens (Thorlabs 
ACT508-200-A-ML). The camera was a Prime BSI Express (Teledyne Photometrics) with 6.5m pixels and 1.8e7 
RMS read noise. No beam expander was used after the detection tube lens. The effective pixel size was 0.130 pm. 
The frame rate was limited by the readout time of the camera. 

The viscosity of this mixture is 7 = 0.0076 Pas, and the density is Pgolvent = 1.18g/mL at T = 22.5°C. The 
expected diffusion coefficient far from the wall is D, = kpT/677nR = 0.0579 pm?/s. The measured average diffusion 
coefficient is 0.0536 ym?/s. Deviations may come from differences in temperature and magnification from nominal 
values and the hydrodynamic interactions with the wall. The density of polystyrene is Ppead = 1.05 g/mL and therefore 
the microspheres are buoyant in this mixture. We observe that the height distribution of the microspheres matches a 
Boltzmann distribution with the expected effective mass m = 40 R3 (Pbead — Psolvent)- 

After reconstructing the sample RI, we identified and localized the microspheres using a custom Python package, 
localize-psf https://github.com/QI2lab/localize-psf. To identify candidate microspheres, we first applied 
a difference-of-Gaussian filter to suppress noise and background, then applied a maximum filter and selected pixels 
where the initial and maximum filtered images had the same values. We keep only points above a threshold. Next, we 
fit an ROI of size 4 x 2 x 2m around each candidate to a 3D Gaussian and only kept spots where the fit parameters 
fell within certain bounds. After localizing the microspheres at all time points, we tracked them using trackpy 
v0.6.1, http: //soft-matter.github.io/trackpy. To mitigate the effect of localization errors, we computed the 
mean-square displacements using every 8th frame. 


F. Preparation of E. coli and microspheres 


A culture of Escherichia coli strain RP437, which is wild type for motility, was inoculated from a single colony and 
incubated overnight at 37°C in Lysogeny Broth (10gL~! tryptone, 5gL~+ yeast extract, 10gL~+ NaCl). Cultures 
were diluted 1:100 in T-broth (10gL~! tryptone, 10¢gL~! NaCl), and incubated at 33°C while shaking at 200 rpm, 
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until reaching an O.Deoo of 0.5. Cells were then centrifuged at 1,200 x g for 7min, washed, and resuspended in 
motility buffer (10mM KzHPO,, 0.1mM EDTA, pH 7.5). The mobility buffer’s RI was 1.333 as measured using a 
refractometer (Kriiss HR 901). The bacterial concentration was estimated to be 1.2 x 10° 1/mL using a hemacytometer 
(Fisher Scientific 0267151B). 

PMS of diameter 0.5 ym (ThermoFisher Fluospheres F8813) of weight/volume 0.02 gmL~' were first sonicated for 
20 min and then diluted by a factor of 3 x 10? with the motility buffer and E. coli mixture. Coverslips were cleaned 
with EtOH. A 120m tall chamber was prepared by placing a secure seal spacer (Electron Microscopy Science Cat 
##70327-20S) on a round #1.5 coverslip of diameter 40mm. ~40 pL of solution were placed in this chamber, and a 
square #1 coverslip with 25mm sides was placed on top to close the chamber. 

To quantify the bacteria swimming and identify the wobble angle, we trained a classifier to distinguish between E. 
coli, tracer particles, and background pixels using LabKit [54]. Next, we labelled individual bacteria and identified 
their main axis using principal component analysis using scikit-image, and tracked the bacteria with trackpy. 
To ensure we tracked only swimming bacteria, we exclude any tracks shorter than 50 frames or which move a total 
distance less than 4pm. To generate the smoothed velocity, we first smoothed the centroid position by computing a 
rolling average over a window of size 15 and then computed the velocity using a second-order difference method [55]. 


G. Preparation of diffusing microspheres in water 


PMS of diameter 1 pm (ThermoFisher Fluospheres F8823) were sonicated for 20min and then 3pL were diluted 
with 97 pL of MQ water. For the shorter chamber, 10 pL of the dilution was pipetted on a #1.5 coverslip of diameter 
40 mm and a square #1 coverslip with 25 mm sides was dropped on top. The chamber was sealed with epoxy (Devcon 
2 Ton Epoxy, GLU-735.90). The height of the fluid chamber was estimated to be ~10 pm from a fluorescence z-stack. 
For the taller chamber, we placed a 120m thick spacer (Electron Microscopy Sciences 70327-20S) on the round 
coverslip, pipetted 40 yL of solution, and closed the chamber with the square coverslip. To achieve a DMD-limited 
frame rate, we cropped the camera chip to 960 x 1040 pixels. We identified and tracked the PMS using the same 
procedure as for the PMS in the water-glycerol mixture. 


SUPPORTING INFORMATION 
$1. FS-ODT MULTIPLEXING VALIDATION WITH MIE THEORY 


We validate our reconstruction approach by comparing with exact results using Mie theory [56]. The Mie theory 
solution to scattering from a spherical dielectric particle expresses the field as a sum of vector spherical harmonics. 
Each spherical harmonic is a product of sinusoids, spherical Bessel functions, and associated Legendre polynomials. 
Although many Python packages calculate the expansion coefficients, few compute the electric fields. We rely on 
miepython [57] to compute the expansion coefficients and developed GPU accelerated Python code to calculate the 
Mie fields. Our code relies on CuPy functions where possible. Since the second order spherical Bessel functions y, 
are not implemented in CuPy, we wrote a CUDA kernel which computes y,, using downward recursion and Miller’s 
algorithm [58]. To avoid the need to rerun this computation for n = 1,...N, our routine provides the full sequence 
Y1,---;Yn- For an image size of 700 x 650 pixels, the CPU implementation runs in ~120s while the GPU runs in ~0.4s. 

Using Mie theory simulations, we validate that our multipled ODT patterns provide considerably more information 
than non-multiplexed patterns and consequently produce superior reconstructions, compared with non-multiplexed 
patterns when the number of patterns is kept fixed. We compare results for 15 images using 10x multiplexing and 
without multiplexing. To generate simulated images, we first compute the Mie electric fields using the expansion 
discussed above for a sphere of diameter 10m with refractive index of 1.59 in media with background refractive 
index 1.515, chosen to match a PMS microspheres in immersion oil. Our routine returns both the plane waves fields 
and the scattered fields. We apply a phase shift to each field, and then sum them to simulate multiplexing. We 
interefere this electric field with a reference plane wave. Then we produced detected fields by applying photon shot 
noise, camera gain and offset, and camera readout noise. We keep the maximum photon number fixed at ~40000 for 
the final images. We reconstruct these simulated images following the same procedure as for the experimental data. 

We find that the multiplexed images (fig. S1A) produce superior reconstructions to the non-multiplexed images 
(fig. S1B). Note that unlike the comparison in fig. 4, the number of images has been kept fixed rather than the 
number of plane waves. Due to the limited number of plane waves used with the non-multiplexed images, significant 
reconstruction artifacts are present in the inferred refractive index. Therefore, for a fixed time budget, we expect that 
multiplexed images provide a superior reconstruction. 
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FIG. $1. Validating multiplexing with Mie theory. A. Maximum intensity projections of refractive index reconstruction 
of a simulated 107m diameter sphere using 10x multiplexing. The axial refractive index reconstruction is slightly distorted 
due to missing cone artifacts and the limited excitation and detection numerical aperture of 1. Scale bar 20pm. B. Maximum 
intensity projections of refractive index reconstruction of a simulated 10 pm diameter sphere without multiplexing. The limited 
angular information produces an inferior reconstruction. 


S2. OPTICAL DIFFRACTION TOMOGRAPHY 


Optical diffraction tomography was performed using a Mach-Zender interferometer, part of a bespoke multimodal 
microscope with quantitative phase imaging (FS-ODT) and multicolor fluorescence superresolution microscopy (struc- 
tured illumination) capabilities. Portions of this microscope are described in previous work [32]. Up to 80mW of 
785 nm light with coherence length of ~50m is generated using a volume-holographic-grating (VHG) stabilized laser 
(Thorlabs FPV785P). This light is divided using a polarizing beam splitter, and the two paths are coupled into sepa- 
rate 1m long polarization-maintaining fibers (Thorlabs PM780-HP). One fiber is used to generate the ODT excitation 
light, and the other is used to generate the reference beam for the off-axis holography. 

The reference beam is collimated with a molded aspheric lens of focal length 13.86 mm (Thorlabs C560TME-B) 
and beam expanded with lenses of focal length 40mm (Thorlabs AC254-040-B-ML) and 300mm (Thorlabs AC508- 
300-AB-ML). It is combined with the reference beam on a d-mirror (Thorlabs BBD1-E03) near the focal plane of the 
40mm lens. 

The excitation light is collimated with a molded aspheric lens of focal length 18.4mm (Thorlabs A280TM-B) and 
beam expanded by lenses with focal length 30mm (Thorlabs AC254-030-AB-ML) and 125mm (Thorlabs LA1986-B- 
ML), resulting in a Gaussian beam with waist ~7mm which is incident on the DMD. 

The ODT patterns are generated from diffraction off a DMD (Texas Instruments DLP6500) with 7.56 ym pitch and 
1920 x 1080 mirrors. We have discussed the details of our DMD geometry elsewhere [32]. Let x point along the long 
axes of the mirror grid and y point along the short axis, both in the plane of the DMD face. Let p = (& + ¥)/V2, 
th = (x — y)/V2, and @ be the normal vector of the DMD surface pointing outwards. The DMD mirrors swivel about 
p the axis and can be in two binary states, either yz = +£12°, which we will call the + and — states. The DMD chip 
is rotated so that the optical table normal points ee p, ensuring that the principle diffraction occurs in the mz 
plane. 

The initial DMD geometry was designed to enable 3 color SIM with excitation wavelength 465 nm, 532nm, and 
635nm. For these three colors to all roughly meet the blaze condition, the DMD face normal z makes an angle of 
04 ~ —21.2° with the optical axis. In the DMD coordinate system, the optical axis points along the sin gm + cos 0g2. 

To achieve high-efficiency ODT, we align the 785nm excitation light to approximately satisfy the blaze condition 
for the (nz,n,) = (—3,3) diffraction order when the mirrors are in the — state. The excitation light is incident on at 
an angle of approximately 6, = 6.64° so that ¥. = siné.m — cos6.%. The light diffracted into (—3,3) travels along 
the 0, = —18.96° along Vv, = sin6,m + cos 62. 

As discussed in the main text, we apply an additional “carrier frequency” to our patterns of frequency f. = 
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+X — ty 1/mirror. This produces additional diffraction about the (—3,3), and we have designed our system such that 
the (—3,3) + (1/4, —-1/4) diffraction order travels along the optical axis and is nearly blazed. The perfectly blazed 
output direction for the — mirrors is sin 6,p+ cos 6,2 for 0, = —17.36°. An aperture blocks all other diffraction orders. 
This specific carrier frequency is chosen to displace the beam as far as possible from orders of the form (n/2,m/2), as 
we see significant diffraction due to the + mirror states. Furthermore, this avoids diffraction along the x and y axes 
coming from a row of DMD mirrors beyond the active chip, which are fixed in the — state. 


The DMD is in a conjugate plane to the back focal plane (Fourier plane) of the excitation and detection objec- 
tives. After light diffracts off of the DMD, it is relayed by a pair of imaging systems, the first using lenses of focal 
length 200 mm (Nikon MXA20696) and 100mm (Thorlabs AC508-100-A-ML), and the second using lenses of 400 mm 
(Thorlabs AC508-400-A-ML) and 300mm (Thorlabs AC508-300-A-ML). After the Nikon tube lens, the NIR light is 
separated from the visible light with a dichroic mirror (Semrock FF750-SDi02-25x36). After the 100mm achromat, 
they are recombined using a second identical dichroic mirror. The two dichroic mirrors are arranged at right angles 
to each other to avoid the differential s— and p—phase shifts from affecting the polarization of the visible light [59]. 
We align the polarization orthogonal to the table surface to avoid similar polarization degradation of the ODT beam 
by the epifluorescence dichroic for the fluorescence modality. 


The ODT light is focused with an oil immersion objective (Olympus UPlanFL N 100x NA 1.3), interacts with the 
sample, and is collected using a water objective (Olympus LUMPLFLN60XW) and a 180mm tube lens (Thorlabs 
AC508-180-AB-ML). The image is then magnified by a factor of 3 using a relay composed of a 100mm (Thorlabs 
AC508-100-B-ML) and 300mm (Thorlabs AC508-300-AB-ML) lenses. The light is imaged onto a Phantom camera 
(VEO-1010L-72G-M) which has 1280 x 960 pixels, pixel size 181m, quantum efficiency ~51% at 785 nm, read-noise 
10.5e— RMS and gain of 0.39 ADU/e~, measured using the approach of [60]. The effective pixel size is 0.1m. The 
maximum frame rate for the full chip is 8420 frame/s, but it can be increased by cropping the chip. 
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FIG. S2. DMD geometry and pattern simulation. A. DMD geometry showing the incident NIR light (red) and the 
main diffraction orders for the + (cyan) and — (red) mirrors and the positive and negative diffraction at the carrier frequencies 
(pale red) B. Magnitude squared of the diffracted electric field versus frequency, |E(f)|?. The diffracted weight is concentrated 
around the positive (red) and negative (magenta) carrier frequency diffraction orders C. Blaze envelope for — mirrors. D. 
Sample DMD pattern with R = 10 mirror show mirrors in the — state (red) and the + state (cyan) E. Pattern in sample plane 
with measured waist wo = 9.31m F. Blaze envelope for the + mirrors. The strength of this diffraction is suppressed by ~10~° 
compared with the — mirrors. 
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A. DMD efficiency 


The expected peak DMD diffraction efficiency into any order (nz, ny) is ~50 %, limited by the reflective of the DMD 
mirrors, the fraction of the chip the mirrors cover, the transmissivity of the DMD window, and diffraction physics. 
Using the carrier frequency, we expect that efficiency into the +1 carrier orders are each ~25% of the power in the 
Oth order. We expect the power in Oth order is ~18% of the power that would be diffracted if the mirrors were all in 
the — state. Taken together, we estimate the DMD diffraction efficiency into the desired order is ~2%. 

The fiber coupling efficiency is ~50%. Additionally, since the fluorescence modality of our microscope operates at 
visible wavelengths, most of the optical coatings are optimized for visible light. Thus, only ~75% of light diffracted 
from the DMD reaches the objectives. The two objective lenses have a combined transmissivity of ~50%. The 
efficiency from the DMD to the camera is thus ~30%, 

As only a small fraction of the DMD mirrors are used to generate a plane wave, this further reduces the efficiency. 
We adjust the magnification between the back focal plane and the DMD so that the pupil radius is approximately the 
same size as the DMD along its narrow dimension. As the magnification factor is M = 0.625, the detection objective 
pupil radius is R, = NA f/0.625 = 4.8mm at the DMD. For a spot pattern of radius R and assuming the laser power 
is distributed uniformly over the pupil, we expect the number of photons per second that strike the camera is 


s0mW rR? 
ope x aR x (0.5 x 0.02 x 0.3) # 2 x 10! photon/s. (S1) 


for R = 10x 7.56m. Here, we divide the terms into incident power, geometric efficiency, and transmission/diffraction 
efficiency. 

For patterns of radius 10 mirrors, our simulations show that the beam waist in the imaging plane is ~10 pm 
(Fig. S2E), and putting this all together, for an imaging time of 100 ps we expect to collect about 


20pm x 755)” 
(20pm x 735)” x 1001s x QE ~ 400 photon (S2) 


N =2x 10!! photon/s x 
4 / x (10 pm)? 


per pixel, where the quantum efficiency of the camera is QE ~ 50% at 785nm. 


B. Pattern fidelity 


Previous ODT approaches using binary DMD patterns have resulted in lower-quality reconstruction than gray- 
scale approaches due to the unwanted additional diffraction orders introduced by the DMD’s square binary pixels. 
Our previous structured illumination microscopy work addressed similar issues [32]. In both cases, these spurious 
diffraction orders arise when using large-scale periodic patterns covering the face of the DMD. In this case, our 
small spot patterns lead to much smaller contributions from these orders, which are mitigated by Fourier broadening 
(Fig. S2B). 


C. DMD non-planarity 


Unlike in many DMD imaging applications, we do not place the DMD orthogonal to the optical axis (Fig. $2A). 
This compromise improves the diffraction efficiency by satisfying the blaze condition. In our geometry, the effect of 
this tilt is minor. 

The tilt introduces a shift in the focus of the plane waves across the DMD face. At the pupil radius, the DMD 
z-shift is at most Rpsin@p ~ 1.73mm. As the z-magnification is M, = M? = 0.39 the shift in the objective back-focal 
plane is 0.67mm. The shift must be compared with the depth-of-focus of the beam in the back focal plane, which we 
estimate using the Rayleigh range. For a plane wave with w, = 10m in the focal plane, the waist in the objective 
BFP is ~37 1m corresponding to a Rayleigh range of tw2/A ~ 5.5mm, which is one order of magnitude larger than 
the focal shift. 

The tilt also introduces deformations in the transformation between the position on the DMD face and the frequency 
of the beam. For example, there is some shear, and a ring pattern on the DMD maps to an oval in frequency space. 
The precise effects can be calculated using the approach of [32] described in section 3 of the supplemental methods. 
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FIG. 53. System stability. A. Exemplary phase stability plot from the first pattern from the 10m bead data. Phases are 
unwrapped B. Exemplary pattern position stability for the same pattern as A. C. Exemplary pattern frequency stability for 
one frequency in the pattern in A. 


D. System stability 


Due to the relatively long beam path (> 2m) used in our multimodal microscope setup, we observe several sources 
of instability in our system and ODT patterns, which we correct for computationally during ODT reconstruction. 
Specifically, we correct for (1) phase drift between the reference arm and the imaging arm, (2) frequency instability 
of the ODT patterns, and (3) position instability of the ODT patterns. 

To correct for (1), we determine the complex factor relating the image electric fields to a single background electric 
field using a least-squares fit. To correct for (2), we determine the location of the Fourier peaks versus time by fitting 
the Fourier transform of the hologram image to a Gaussian in the vicinity of each peak. To correct for (3), we register 
images using the Fourier transform of the absolute value of the electric field. Empirically, we find that using the 
absolute value of the electric field is superior to using either the Fourier transform of the intensity or the Fourier 
transform of the electric field. The improvement is likely due to the rejection of background fluctuations that do not 
involve the interference pattern and taking the absolute value removes the need to consider phase fluctuations of the 
electric field. 

The stability of our system over ~4s is shown in Fig. $3. Over this short time span the relative phase between the 
imaging and reference beam drifts by ~ 27 radians, the spatial position drifts by <50 nm, and the beam frequency drifts 
by <4x 107? pm’. Note that the Fourier space pixel sizes are df, = 7.81 x 107? ym™ and df, = 10.42 x 107? pm™?, 
so the frequency drift is sub-pixel. 


$3. ITERATIVE RECONSTRUCTION WITH PROXIMAL GRADIENT METHODS 


As FS-ODT multiplexing introduces a more challenging computational image reconstruction problem requiring 
new approaches we briefly discuss some commonly used reconstruction approaches. In 1969, Wolf introduced optical 
diffraction tomography using a Born approximation formulation [8]. In this formulation, the forward model describing 
the electric field after interacting with a given RI distribution is linear. This linearity makes solving the inverse 
problem, i.e. inferring the RI based on the observed electric fields, particularly simple. Devaney realized that the 
Rytov approximation is better suited for biological sample [9], and various technical improvements have expanded 
the range of validity and quality of reconstructions [61-68]. However, imaging thicker and higher-contrast samples 
inevitably introduces multiple scattering, which the linear Born and Rytov approximations do not include. To 
address multiple scattering, a variety of multi-slice reconstruction approaches have been developed that rely on 
the beam-propagation model (BPM) [69-73]. However, the BPM entails the paraxial approximation, motivating 
the development of more accurate forward models, including the split-step non-paraxial (SSNP) model [36], more 
sophisticated Born approximation approaches [74-76], and HyPM [77]. Most approaches consider forward scattering 


19 


only, which is implicit in the layer-by-layer multi-slice forward models. Other approaches account for backscattering 
using the Lippman-Schwinger equation [78-80], but these are generally more computationally expensive than multi- 
slice models. Machine learning approaches are increasingly employed to accelerate and denoise RI reconstruction [44]. 
In this work, we primarily rely on the BPM or the SSNP combined with an initial guess generated using a demultiplexed 
low-resolution Rytov approximation approach to infer RI information from FS-ODT data. 

We treat refractive index reconstruction as a regularized minimization problem and solve it using the fast iterative 
shrinkage-thresholding algorithm (FISTA) [35]. As usual in FISTA, we attempt to minimize a function which is 
the sum of two terms, a loss function and a regularization function L(x) + g(x), where CL is differentiable and has a 
Lipschitz-continuous gradient with Lipschitz constant L, and g is convex. We iteratively update the proposed solution, 
Zz, starting from t = 0, go = 1, the step-size y, and a starting guess xr 


Yt = prox, [vr-1 — YVL(a1-1)] 


1+,/1+4q?_, 


t= 5} 
G@-1—1 
Tt+1 = ye + <——— (He — yea). 
t 


In the last step the convergence is accelerated by including a momentum term and the structure of q@ is chosen to 
change the convergence from O(1/t) to O(1/t?). 
The proximal operator for g, 


1 
prox, (2) = argmin, {g(e) + 5 Je 2b, (S3) 
27 


determines a new object which is near to the initial value but better satisfies the regularization. When an explicit 
form of or fast algorithm for computing the proximal operator is known, this process is efficient, as in the case of 
total variation [81], ¢:, or 2 norms. In this work we choose a regularization function which enforces smoothness and 
sparsity 


g(x) = TevTV (a) + Te, || II1. (S4) 


The total loss, £(a,) + g(#z), is non-increasing when the step-size y € [0,2/LZ]. If the Lipschitz constant is not 
known, a line-search strategy can determine the step-size at each iteration by reducing an initial step-size until the 
Lipschitz condition is locally satisfied (algorithm 1). This requires additional calculation of L(y), £(x), and VL(y). 


Data: Proposed object x 
Data: Initial step-size yo 
Data: Step-size adjuster a < 1 
Y= Yo3 
y © prox [x — yVL(x)]; 
while L(y) > L(x) +VL()-(y— 2) + lly — x||3 do 

yo ay; 

y © prox [x — yVL(x)]; 
end 

Algorithm 1: Line-search algorithm for setting the step-size 


In ODT, we measure a sequence of electric fields B®, where i = 1,...,.M indexes the incident angles, derived from 
off-axis holography and define our loss function by 


L= flLet+(l—f)£r (S5) 
x(n) = sz ow — wf (50 
Lil 1 AS gO] _[y@al]? 

i(n) = sa 2 2||¥ Ie | cD 
== 


where f € [0,1] describes the relative strength of the phase-sensitive (C~) and phase-insensitive (£;) components of 
the loss function, &)(n) is the forward model describing the predicted electric field as a function of n, || - || is the 02 
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norm, and N is the number of pixels in a single electric field. Note that many different loss functions can be defined 
for the phase-insensitive portion, and in fact the most natural one might seem to involve |W|?, however previous work 
has demonstrated better convergence using eq. S7 [82]. 

To perform FISTA, we require the gradient of the loss function with respect to the refractive index. The loss 
function, £ : RN x RN — R, can be interpreted either using the real and imaginary parts of the refractive index 
as components of a vector, or as complex numbers. Since £ is not holomorphic, if we regard the domain as CN we 
must compute the gradient using Wirtinger derivatives [34], defined by 0, = $ (0, — 10.) and 0,» = $ (0, + 102), 
where z’ and 2” are the real and imaginary parts of z respectively. In this formalism, the gradient used in FISTA is 
V = 20,». Computing this derivative, which only acts directly on the forward model, we find (suppressing the beam 
angle index) 


sme fe . 
Atha x [w = w| (S9) 
Awe walbd = Jel] o ww (S10) 


where b and / index the xy- and z-position of the voxels respectively and we have used the chain rule 0,» (£0 V) = 
On«L (On V)* + OnL (On«V) and the fact On«V = 0. 

Due to the FS-ODT pattern generation strategy, our illumination beams typically do not cover the entire field of 
view and exhibit decreasing intensity near the edges. In regions with less intensity, the influence of the loss function 
decreases and the effect of the regularization increases. Therefore, to achieve a higher quality reconstruction over a 
larger field of view it is helpful to normalize the loss function by the electric field. We optionally replace the loss 
function with 


eee. =e 
(S11) 
te= Lae ws sey? 
xo atl -bell 
ey (S12) 


i=1 j=1 ws] + aq? 


where a is a regularization parameter that prevents division by small numbers where the electric field is noise domi- 
nated. The denominator does not depend on the forward model, so the loss function gradients need only be divided 
elementwise by the denominator to account for this change. 


A. Linear scattering models 


As usual we suppose light interacts with a spatially varying refractive index according to the scalar Helmholtz 
equation and work with phasors carrying exp(iwt) time dependence, 


[V? + k?n?(r)] E(r) = 0. (S13) 


where k = 20/Ano. 
Additionally we define the scattering potential V 


on\2 
V(r) =- (=) [n2(r) — n?] (S14) 

We suppose that our sample is illuminated by a sequence of plane waves and the ith plane wave has frequency f 
E%°) = exp [-iane : r| (S15) 


where 27 [£| =k. 
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In the Born approximation, valid when the cumulative phase shift of the beam is < 7/4 [83], the 2D Fourier 
transform of the scattered electric field gives the scattering potential along a spherical cap in 3D Fourier space 


F(i,s) = a ; @) ¢ — ¢@ ¢ — pO) (Born 
E (fe, fy) = Dx Ofte) re ity 5 itz ce ) (Bo ) (S16) 


The Rytov approximation is an alternate approach which is usually more accurate for biological samples. In this 
approximation, the scattering potential is related to the Rytov phase w(r), 


BE) (p) = E,(r) Cae be 1) (S17) 
. (@ 
~(r) = log Pare +74 unwrap {angle Eg (r)] — angle [ ete) (r)] \ (S18) 
VO fe I Fy 49 y= : VGe Je sigh) set) (Rysey) (S19) 


Dix On filfe fy) 


The Rytov approximation is valid when n?(r) — n2 >> |Vy(xr)|? k? [9]. 
In this case, it is more convenient to work with the scattering potential than the refractive index, and the loss 
function (eq. 56) and its gradient are 


where &) is the scattered field or the Rytov phase depending on the approximation, F™ is the forward model linear 
operator for the 7th angle which connects the 3D Fourier transform of the scattering potential, V, to ¥, The second 
factor of 1/N converts this to the loss in real-space accounting for the Fourier transform. Since the forward model 
is linear, a Lipschitz constant for the loss function is proportional to the largest eigenvalue of F'F which can be 
computed using a singular value decomposition or the power iteration algorithm. 


B. Multi-slice models 


Suppose that we can rewrite the Helmholtz equation (eq. $13) by reparameterizing the electric field as V and 
rewriting the differential operator as the sum of two terms, where A describes the effect of the background refractive 
index n, and B captures the effect of spatially varying refractive index perturbations, 


0.V = (A+ B)W. (S20) 
Formally the solution is a path-ordered exponential, but for small enough 62, 


W(dz) = exp [(A + B)dz] U(0 (S21) 
= exp [Adz] exp [Bdz] (0). (S22) 


For convenience we define P = exp [A6z] and Q = exp[Béz]. Corrections to eq. $22 are given by the Baker-Campbell- 
Hausdorff formula, and the first correction is proportional to the commutator [A, B]. By construction, this vanishes 
when n = no, and we expect it and higher order commutators to be small when the refractive index perturbation is 
small. 

We can propagate an initial field through a volume by discretizing it into layers of thickness 6z and iteratively 
applying eq. $22 


U = F(PQ)*-9.. (PQ)VOGO (S23) 


Here F describes model operations beyond the final refractive index plane, ¥(”) is the intermediate field before the 
mth voxel, and W is the detected field. 
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For models considered here, Q is local in the sense that the only dependence on n(z;) is in Q™, and thus 


ov dQ 

a _ (wen pw v) S24 
Onv(Z1) ( Onv(Z1) a ( ) 
WED = F(PQ)E-Y_ (PQ), (S25) 


where the index a in the first equation indicates the ath component of the vector. Typically the structure of OQ allow 
us to simplify this expression. For example, in the BPM Q is diagonal due to the fact it does not mix the field at 
different r positions. For the SSNP, Q is local but mixes the derivative and field at the same position. 


C. Beam-propagation model (BPM) 


In the BPM, the Helmholtz equation is put in the form of eq. $20 by first making the paraxial approximation [71]. 
Here we take V = EF and the model is defined by 


Pv =F! {exp [ise [2 — k2 — k2] x Fy {WV} (he, ky) } (S26) 


now aces ita 6e Mena) — Mo ay (827) 
7) 


PUSS? {H(ke, ky) X exp lider [2 — 2 — ke] x Fi {WV} (ke, hy) } (S28) 


where H is the coherent transfer function, dz, is the final distance the beam propagates, and k, = 27/X. 77 is the 
obliquity factor which is taken to be 1 in most cases, but taking 7 = 1/cos@ improves accuracy [34]. 
Following eqs. $8 and $24 the derivatives are 


@. & 7 
OQ ed = Woh exp med x Ocd x Obd (S29) 
Onv(21) n 
OWe tk d 2 41) (4D) 
= ot? U 
ay oe ($30) 
ikydz t * 
nee (+1) : (+1) 
Vince ; (ow ) Ow c| © (w ) (S31) 


D. Split-step non-paraxial model (SSNP) 


This model is extensively discussed elsewhere [34, 36] and we briefly discuss it here for convenience. To rewrite the 
Helmholtz equation in the form of eq. $20 we must work with a vector of the electric field and its first derivative. The 
physical intuition behind this is as follows. Suppose we know the electric field in a single plane and wish to know its 
value everywhere. We can decompose the field into Fourier modes, but we cannot distinguish forward and backwards 
travelling plane waves at the same lateral spatial frequency. We need additional information to untangle these two 
contributions: e.g. the magnetic field or the axial derivative of the electric field 0,E. Therefore the propagation 
operator must act on both the field and its derivative. 

In this case the forward model is 


E 
oh (ain) (S32) 
= et cos kz6z p sink.6z 
nae Gs GHses Scoskeo en ke ee (S33) 
Oy — 1 0 
iba (1: [n2 — n(x, y, 2)] 62 i) x (S34) 
FU=FO'{(3 376) H (kes ky) x P(dz4)¥}. ($35) 


When we write Q as a matrix following eq. S24 we combine the structure of the r, index and the field/derivative 
index. Let i and j index the position r, and field/derivative respectively and define the composite index c = 2i + j. 
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Following eqs. $8, $24, and $25, the derivatives are 


wee = —2k? Oz nx b¢,2b41 64,26 (S36) 
OEa _ _ opp (41) gl) 

Ono(z1) ake Jz Mm X Yo ones Yop (S37) 
o£ 2 (1+1) (l)\* « 

Ony(z1) ake bz d, as sen in, a (YS) ny (21) (S38) 
Vnl = —2k? 6z [(e cy) 'au-£| © n*(z) © (no) : (S39) 


where Y(+1) = W+) P© and in the last line we use V2, = E, and define Y, a modified version of Y which retains 
only the derivative parts of the second index. 
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